Read in the data

Due to the inconsistent column naming covention, we manually convert the cx column to cX in an effort to reduce obfuscation.

# Read in the data
data <- read.csv('Greatest_Aussie_Groceries_sales_data.csv', header=TRUE, sep=",")
# Change column name of cx to match style of capital X and Y
colnames(data)[colnames(data)=="cx"] <- "cX"

Reformatting Data

# Organic X
dataR <- data %>% filter(class=="organic") %>% select(STORE, WEEK, Hval_150, pX, cX, oz_X, deal_X)
dataR$brand <- "X"
dataR$class <- "organic"
colnames(dataR) <- c("STORE", "WEEK", "Hval", "p", "c", "oz", "deal", "brand", "class")
data.ref <- dataR
# Organic Y
dataR <- data %>% filter(class=="organic") %>% select(STORE, WEEK, Hval_150, pY, cY, oz_Y, deal_Y)
dataR$brand <- "Y"
dataR$class <- "organic"
colnames(dataR) <- c("STORE", "WEEK", "Hval", "p", "c", "oz", "deal", "brand", "class")
data.ref <- rbind(data.ref, dataR)
# NonOrganic X
dataR <- data %>% filter(class=="nonorganic") %>% select(STORE, WEEK, Hval_150, pX, cX, oz_X, deal_X)
dataR$brand <- "X"
dataR$class <- "nonorganic"
colnames(dataR) <- c("STORE", "WEEK", "Hval", "p", "c", "oz", "deal", "brand", "class")
data.ref <- rbind(data.ref, dataR)
# NonOrganic Y
dataR <- data %>% filter(class=="nonorganic") %>% select(STORE, WEEK, Hval_150, pY, cY, oz_Y, deal_Y)
dataR$brand <- "Y"
dataR$class <- "nonorganic"
colnames(dataR) <- c("STORE", "WEEK", "Hval", "p", "c", "oz", "deal", "brand", "class")
data.ref <- rbind(data.ref, dataR)

Mutate Data

We mutate the data to include columns for deal_feat (an indictator for both deal and features), revenue, and profit.

# Append a deal_feat column for X and Y
data <- mutate(data, deal_feat_Y = deal_Y*10 + feat_Y, deal_feat_X = deal_X*10 + feat_X)
# Append a revenue column for X and Y 
data <- mutate(data, rev_X = oz_X * pX, rev_Y = oz_Y * pY)
# Append a profit column for X and Y 
data <- mutate(data, profit_X = rev_X - cX * pX, profit_Y = rev_Y - cY * pY)

Summary Plots

Lets start by just plotting the profit over the 52 weeks for each store.

plotStore <- function(store_ID) {
# The palette with black:
cbbPalette <- c("red4", "orange1", "blue4", "skyblue1", "green4", "chartreuse1", "darkorchid4", "mediumorchid1")
# Pull data into temp results dataframe
results <- data.frame(WEEK=c(1:52))
results[,c("profit_X_org","profit_Y_org")] <- data %>% filter(.,STORE==store_ID, class=="organic") %>% select(.,profit_X, profit_Y)
results[,c("profit_X_non","profit_Y_non")] <- data %>% filter(.,STORE==store_ID, class=="nonorganic") %>% select(.,profit_X, profit_Y)
results[,c("profit_X_org_col","profit_Y_org_col")] <- data %>% filter(.,STORE==store_ID, class=="organic") %>% select(.,deal_X, deal_Y)
results[,c("profit_X_non_col","profit_Y_non_col")] <- data %>% filter(.,STORE==store_ID, class=="nonorganic") %>% select(.,deal_X, deal_Y)
# Assign legend name to categorical data
results$profit_X_org_col[results$profit_X_org_col == 0] <- "X organic w/o deal"
results$profit_X_org_col[results$profit_X_org_col == 1] <- "X organic w/ deal"
results$profit_Y_org_col[results$profit_Y_org_col == 0] <- "Y organic w/o deal"
results$profit_Y_org_col[results$profit_Y_org_col == 1] <- "Y organic w deal"
results$profit_X_non_col[results$profit_X_non_col == 0] <- "X nonorganic w/o deal"
results$profit_X_non_col[results$profit_X_non_col == 1] <- "X nonorganic w deal"
results$profit_Y_non_col[results$profit_Y_non_col == 0] <- "Y nonorganic w/o deal"
results$profit_Y_non_col[results$profit_Y_non_col == 1] <- "Y nonorganic w deal"
# Plot results
ggplot(results, aes(x=WEEK)) + 
  geom_point(aes(y=profit_X_org, colour=profit_X_org_col)) +
  geom_point(aes(y=profit_Y_org, colour=profit_Y_org_col)) + 
  geom_point(aes(y=profit_X_non, colour=profit_X_non_col)) +
  geom_point(aes(y=profit_Y_non, colour=profit_Y_non_col)) + 
  scale_colour_manual(values=cbbPalette) +
  labs(x = "Week", y = "Profit", title = paste("Profit for Store",store_ID)) 
}

# Plot the data
plotStore(1)

plotStore(2)

plotStore(3)

plotStore(4)

plotStore(5)

plotStore(6)

plotStore(7)

Box plots of profit for the different products

This is just a simple box plot analysis

# Function to retreve data from dataframe and composite into factor
retrieveData <- function(data, deal, CLASS, xy) {
  names <- c("STORE", "PROFIT")
  if (xy == "x")
    temp <- data %>% filter(deal_X==deal, class==CLASS) %>% select(STORE, profit_X)
   else
    temp <- data %>% filter(deal_Y==deal, class==CLASS) %>% select(STORE, profit_Y)
  colnames(temp) <- names
  name <- if (xy == "x") "X" else "Y"
  name <- if (CLASS == "organic") paste(name,"O",sep="") else paste(name,"I",sep="") 
  name <- if (deal == 1) paste(name,"deal",sep="_") else paste(name,"no_deal",sep="_") 
  return(data.frame(type=rep(name,nrow(temp)),temp))
}

boxplot_data <- retrieveData(data, deal=0, CLASS="organic", xy="x")
boxplot_data <- rbind(boxplot_data,retrieveData(data, deal=0, CLASS="nonorganic", xy="x"))
boxplot_data <- rbind(boxplot_data,retrieveData(data, deal=0, CLASS="organic", xy="y"))
boxplot_data <- rbind(boxplot_data,retrieveData(data, deal=0, CLASS="nonorganic", xy="y"))
boxplot_data <- rbind(boxplot_data,retrieveData(data, deal=1, CLASS="organic", xy="x"))
boxplot_data <- rbind(boxplot_data,retrieveData(data, deal=1, CLASS="nonorganic", xy="x"))
boxplot_data <- rbind(boxplot_data,retrieveData(data, deal=1, CLASS="organic", xy="y"))
boxplot_data <- rbind(boxplot_data,retrieveData(data, deal=1, CLASS="nonorganic", xy="y"))

ggplot(boxplot_data, aes(x=type, y=PROFIT)) + geom_boxplot() + labs(title="Boxplots of profit for different products", x="Product", y = "Profit")

Given the massive spread between no deal and deal data, let’s take a log10() scale transform of the y-axis

boxplot_data <- mutate(boxplot_data, log10PROFIT=log10(PROFIT))
ggplot(boxplot_data, aes(x=type, y=log10PROFIT)) + geom_boxplot() + labs(title="Boxplots of log10(profit) for different products", x="Product", y = "log10(Profit)")

It is clear from these boxplots that the median values for X and Y products are approximately equal. The spread of profit for Y is much larger than X. Profit increases dramatically when a deal is going on.

Just as an additional spam of figures, lets look at the box plots for each product seperated by store (again with a log10 scaling)

boxplot_data %>% filter(type=="XO_no_deal") %>% ggplot(., aes(x=STORE, y=PROFIT)) + geom_boxplot(aes(group = cut_width(STORE, 1))) + labs(title=paste("Profit of", "XO_no_deal" ,"for different stores"), x="Store", y = "Profit")

boxplot_data %>% filter(type=="XI_no_deal") %>% ggplot(., aes(x=STORE, y=PROFIT)) + geom_boxplot(aes(group = cut_width(STORE, 1))) + labs(title=paste("Profit of", "XI_no_deal" ,"for different stores"), x="Store", y = "Profit")

boxplot_data %>% filter(type=="YO_no_deal") %>% ggplot(., aes(x=STORE, y=PROFIT)) + geom_boxplot(aes(group = cut_width(STORE, 1))) + labs(title=paste("Profit of", "YO_no_deal" ,"for different stores"), x="Store", y = "Profit")

boxplot_data %>% filter(type=="YI_no_deal") %>% ggplot(., aes(x=STORE, y=PROFIT)) + geom_boxplot(aes(group = cut_width(STORE, 1))) + labs(title=paste("Profit of", "YI_no_deal" ,"for different stores"), x="Store", y = "Profit")

boxplot_data %>% filter(type=="XO_deal") %>% ggplot(., aes(x=STORE, y=PROFIT)) + geom_boxplot(aes(group = cut_width(STORE, 1))) + labs(title=paste("Profit of", "XO_deal" ,"for different stores"), x="Store", y = "Profit")

boxplot_data %>% filter(type=="XI_deal") %>% ggplot(., aes(x=STORE, y=PROFIT)) + geom_boxplot(aes(group = cut_width(STORE, 1))) + labs(title=paste("Profit of", "XI_deal" ,"for different stores"), x="Store", y = "Profit")

boxplot_data %>% filter(type=="YO_deal") %>% ggplot(., aes(x=STORE, y=PROFIT)) + geom_boxplot(aes(group = cut_width(STORE, 1))) + labs(title=paste("Profit of", "YO_deal" ,"for different stores"), x="Store", y = "Profit")

boxplot_data %>% filter(type=="YI_deal") %>% ggplot(., aes(x=STORE, y=PROFIT)) + geom_boxplot(aes(group = cut_width(STORE, 1))) + labs(title=paste("Profit of", "YI_deal" ,"for different stores"), x="Store", y = "Profit")

Regression

# Function that pulls out the data
retrieveDataLM <- function(data, CLASS, xy) {
  names <- c("STORE", "oz", "p", "deal")
  if (xy == "x")
    temp <- data %>% filter(.,class==CLASS) %>% select(.,STORE, oz_X, pX, deal_X)
  else
    temp <- data %>% filter(class==CLASS) %>% select(STORE, oz_Y, pY, deal_Y)
  colnames(temp) <- names
  name <- if (xy == "x") "X" else "Y"
  name <- if (CLASS == "organic") paste(name,"O",sep="") else paste(name,"I",sep="") 
  return(data.frame(type=rep(name,nrow(temp)),temp))
}

# Function that plots the loglog and actual curves
filteredLM <- function(data, TYPE, log=FALSE, include_Deal=FALSE) {
  response <- data %>% filter(type==TYPE) %>% select(oz) %>% {if (log) log(.) else (.)} %>% as.matrix()
  elasticity <- data %>% filter(type==TYPE) %>% select(p) %>% {if (log) log(.) else (.)} %>% as.matrix()
  deal <- data %>% filter(type==TYPE) %>% select(deal) %>% {if (include_Deal) (.) else (.)*0} %>% as.matrix()
  return(lm(response~elasticity+deal))
}

# Function to plot Prediction with actual data
plotStatsLogLog <- function(lm, type="loglog", title="") {
  results <- data.frame(p=lm$model$elasticity, oz=lm$model$response, resid=resid(lm))
  results <- mutate(results, loglogSol=lm$coefficients[1] + lm$coefficients[2]*p + if (is.na(lm$coefficients[3])) 0 else lm$coefficients[3]*lm$model$deal)
  results$sol <- exp(results$loglogSol)
  if (type == "loglog")
    ggplot(results) + geom_point(aes(x=p, y=oz)) + geom_line(aes(x=p, y=loglogSol)) + labs(x="log(price)", y="log(oz)", title=title) + coord_flip()
  else
    ggplot(results) + geom_point(aes(x=exp(p), y=exp(oz))) + geom_line(aes(x=exp(p), y=sol))  + labs(x="price", y="oz", title=title) + coord_flip()
}

# Function that generates the inear model
lm_data <- retrieveDataLM(data, CLASS="organic", xy="x")
lm_data <- rbind(lm_data,retrieveDataLM(data, CLASS="nonorganic", xy="x"))
lm_data <- rbind(lm_data,retrieveDataLM(data, CLASS="organic", xy="y"))
lm_data <- rbind(lm_data,retrieveDataLM(data, CLASS="nonorganic", xy="y"))

Log Regression for Organic X without deal dummy variable

Linear regression for \(\log(oz_{x-org}) = \log(k) + r \log(p_{x-org})\)

# Linear Model Summary for Organic X product without including the deal dummy variable
XO_lm_without_deal <- filteredLM(lm_data, TYPE="XO", log=TRUE, include_Deal=FALSE)
summary(XO_lm_without_deal)
## 
## Call:
## lm(formula = response ~ elasticity + deal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.65001 -0.23331 -0.06901  0.14514  1.17942 
## 
## Coefficients: (1 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -16.8432     0.5183  -32.49   <2e-16 ***
## elasticity   -7.0398     0.1445  -48.73   <2e-16 ***
## deal              NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3749 on 362 degrees of freedom
## Multiple R-squared:  0.8677, Adjusted R-squared:  0.8674 
## F-statistic:  2375 on 1 and 362 DF,  p-value: < 2.2e-16
plotStatsLogLog(XO_lm_without_deal, "loglog", "loglog PED for Organic X without deal")

plotStatsLogLog(XO_lm_without_deal, "normal", "PED for Organic X without deal")

Log Regression for Organic X with deal dummy variable

Linear regression for \(\log(oz_{x-org}) = \log(k) + r \log(p_{x-org}) + deal_{x-org}\)

# Linear Model Summary for Organic X product including deal dummy variable
XO_lm_with_deal <- filteredLM(lm_data, TYPE="XO", log=TRUE, include_Deal=TRUE)
summary(XO_lm_with_deal)
## 
## Call:
## lm(formula = response ~ elasticity + deal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.48496 -0.12887 -0.03184  0.08066  1.21545 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.16811    0.71250  -4.446 1.16e-05 ***
## elasticity  -3.12964    0.20256 -15.450  < 2e-16 ***
## deal         1.39566    0.06387  21.851  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2463 on 361 degrees of freedom
## Multiple R-squared:  0.9431, Adjusted R-squared:  0.9427 
## F-statistic:  2989 on 2 and 361 DF,  p-value: < 2.2e-16
plotStatsLogLog(XO_lm_with_deal, "loglog", "loglog PED for Organic X with deal")

plotStatsLogLog(XO_lm_with_deal, "normal", "PED for Organic X with deal")

Log Regression for Organic Y without deal dummy variable

Linear regression for \(\log(oz_{y-org}) = \log(k) + r \log(p_{y-org})\)

# Linear Model Summary for Organic Y product without including the deal dummy variable
YO_lm_without_deal <- filteredLM(lm_data, TYPE="YO", log=TRUE, include_Deal=FALSE)
summary(YO_lm_without_deal)
## 
## Call:
## lm(formula = response ~ elasticity + deal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.85175 -0.65452  0.03862  0.62843  2.76432 
## 
## Coefficients: (1 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -15.7803     1.2829  -12.30   <2e-16 ***
## elasticity   -6.7308     0.3595  -18.72   <2e-16 ***
## deal              NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9754 on 362 degrees of freedom
## Multiple R-squared:  0.4919, Adjusted R-squared:  0.4905 
## F-statistic: 350.5 on 1 and 362 DF,  p-value: < 2.2e-16
plotStatsLogLog(YO_lm_without_deal, "loglog", "loglog PED for Organic Y without deal")

plotStatsLogLog(YO_lm_without_deal, "normal", "PED for Organic Y without deal")

Log Regression for Organic Y with deal dummy variable

Linear regression for \(\log(oz_{y-org}) = \log(k) + r \log(p_{y-org}) + deal_{y-org}\)

# Linear Model Summary for Organic Y product including deal dummy variable
YO_lm_with_deal <- filteredLM(lm_data, TYPE="YO", log=TRUE, include_Deal=TRUE)
summary(YO_lm_with_deal)
## 
## Call:
## lm(formula = response ~ elasticity + deal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.79111 -0.65105  0.03551  0.60161  2.80730 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -5.9350     2.3642  -2.510   0.0125 *  
## elasticity   -3.9058     0.6740  -5.795 1.49e-08 ***
## deal          1.1975     0.2445   4.897 1.47e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9458 on 361 degrees of freedom
## Multiple R-squared:  0.5236, Adjusted R-squared:  0.5209 
## F-statistic: 198.4 on 2 and 361 DF,  p-value: < 2.2e-16
plotStatsLogLog(YO_lm_with_deal, "loglog", "loglog PED for Organic Y with deal")

plotStatsLogLog(YO_lm_with_deal, "normal", "PED for Organic Y with deal")

Log Regression for Nonorganic X without deal dummy variable

Linear regression for \(\log(oz_{x-nonorg}) = \log(k) + r \log(p_{x-nonorg})\)

# Linear Model Summary for Nonorganic X product without including the deal dummy variable
XI_lm_without_deal <- filteredLM(lm_data, TYPE="XI", log=TRUE, include_Deal=FALSE)
summary(XI_lm_without_deal)
## 
## Call:
## lm(formula = response ~ elasticity + deal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.62446 -0.23014 -0.09199  0.08494  1.27804 
## 
## Coefficients: (1 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -15.8704     0.5801  -27.36   <2e-16 ***
## elasticity   -6.7511     0.1630  -41.43   <2e-16 ***
## deal              NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3569 on 362 degrees of freedom
## Multiple R-squared:  0.8258, Adjusted R-squared:  0.8253 
## F-statistic:  1716 on 1 and 362 DF,  p-value: < 2.2e-16
plotStatsLogLog(XI_lm_without_deal, "loglog", "loglog PED for Nonorganic X without deal")

plotStatsLogLog(XI_lm_without_deal, "normal", "PED for Nonrganic X without deal")

Log Regression for Nonorganic X with deal dummy variable

Linear regression for \(\log(oz_{x-nonorg}) = \log(k) + r \log(p_{x-nonorg}) + deal_{x-nonorg}\)

# Linear Model Summary for Nonorganic X product including deal dummy variable
XI_lm_with_deal <- filteredLM(lm_data, TYPE="XI", log=TRUE, include_Deal=TRUE)
summary(XI_lm_with_deal)
## 
## Call:
## lm(formula = response ~ elasticity + deal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.52573 -0.10589 -0.00324  0.06946  1.30459 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.27234    0.60919  -5.372  1.4e-07 ***
## elasticity  -3.15079    0.17312 -18.200  < 2e-16 ***
## deal         1.40060    0.05548  25.243  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2149 on 361 degrees of freedom
## Multiple R-squared:  0.937,  Adjusted R-squared:  0.9367 
## F-statistic:  2685 on 2 and 361 DF,  p-value: < 2.2e-16
plotStatsLogLog(XI_lm_with_deal, "loglog", "loglog PED for Nonorganic X with deal")

plotStatsLogLog(XI_lm_with_deal, "normal", "PED for Nonorganic X with deal")

Log Regression for Nonorganic Y without deal dummy variable

Linear regression for \(\log(oz_{y-nonorg}) = \log(k) + r \log(p_{y-nonorg})\)

# Linear Model Summary for Nonorganic Y product without including the deal dummy variable
YI_lm_without_deal <- filteredLM(lm_data, TYPE="YI", log=TRUE, include_Deal=FALSE)
summary(YI_lm_without_deal)
## 
## Call:
## lm(formula = response ~ elasticity + deal)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5546 -0.6087  0.0357  0.6109  3.4261 
## 
## Coefficients: (1 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -15.4133     1.2595  -12.24   <2e-16 ***
## elasticity   -6.6570     0.3529  -18.86   <2e-16 ***
## deal              NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9586 on 362 degrees of freedom
## Multiple R-squared:  0.4957, Adjusted R-squared:  0.4943 
## F-statistic: 355.8 on 1 and 362 DF,  p-value: < 2.2e-16
plotStatsLogLog(YI_lm_without_deal, "loglog", "loglog PED for Nonorganic Y without deal")

plotStatsLogLog(YI_lm_without_deal, "normal", "PED for Nonorganic Y without deal")

Log Regression for Nonorganic Y with deal dummy variable

Linear regression for \(\log(oz_{y-nonorg}) = \log(k) + r \log(p_{y-nonorg}) + deal_{y-nonorg}\)

# Linear Model Summary for Nonorganic Y product including deal dummy variable
YI_lm_with_deal <- filteredLM(lm_data, TYPE="YI", log=TRUE, include_Deal=TRUE)
summary(YI_lm_with_deal)
## 
## Call:
## lm(formula = response ~ elasticity + deal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.46784 -0.56287  0.02177  0.56931  3.02660 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.3106     2.3526  -0.982    0.327    
## elasticity   -2.8956     0.6712  -4.314 2.07e-05 ***
## deal          1.5425     0.2386   6.464 3.31e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9087 on 361 degrees of freedom
## Multiple R-squared:  0.548,  Adjusted R-squared:  0.5455 
## F-statistic: 218.8 on 2 and 361 DF,  p-value: < 2.2e-16
plotStatsLogLog(YI_lm_with_deal, "loglog", "loglog PED for Nonorganic Y with deal")

plotStatsLogLog(YI_lm_with_deal, "normal", "PED for Nonorganic Y with deal")

Cross Elasticity of demand

# Function that pulls out the data
retrieveDataCPED <- function(data, Q_class, Q_brand) {
tempQ <- data %>% filter(., class==Q_class) %>% select(., STORE, matches(paste("oz",Q_brand,sep="_")))
colnames(tempQ) <- c("STORE", "oz")
tempP_Org <- data %>% filter(., class=="organic") %>% select(., pX, pY)
colnames(tempP_Org) <- c("p_X_org", "p_Y_org")
tempP_NonOrg <- data %>% filter(., class=="nonorganic") %>% select(., pX, pY)
colnames(tempP_NonOrg) <- c("p_X_nonorg", "p_Y_nonorg")
return(data.frame(tempQ, tempP_Org, tempP_NonOrg))
}

# Function that plots the loglog and actual curves
crossPEDLinearReg <- function(data, TYPES, log=TRUE) {
  Q <- data %>% select(oz) %>% {if (log) log(.) else (.)} %>% as.matrix()
  pXOrg <- data %>% select(p_X_org) %>% {if (log) log(.) else (.)} %>% {if ("xOrg" %in% TYPES) (.) else (.)*0} %>% as.matrix()
  pYOrg <- data %>% select(p_Y_org) %>% {if (log) log(.) else (.)} %>% {if ("yOrg" %in% TYPES) (.) else (.)*0}  %>% as.matrix()
  pXNonOrg <- data %>% select(p_X_nonorg) %>% {if (log) log(.) else (.)} %>% {if ("xNonOrg" %in% TYPES) (.) else (.)*0}  %>% as.matrix()
  pYNonOrg <- data %>% select(p_Y_nonorg) %>% {if (log) log(.) else (.)} %>% {if ("yNonOrg" %in% TYPES) (.) else (.)*0}  %>% as.matrix()
  return(lm(Q~pXOrg+pYOrg+pXNonOrg+pYNonOrg))
}


crossXOrg <- retrieveDataCPED(data, "organic", "X")
crossYOrg <- retrieveDataCPED(data, "organic", "Y")
crossXNonOrg <- retrieveDataCPED(data, "nonorganic", "X")
crossYNonOrg <- retrieveDataCPED(data, "nonorganic", "Y")

lmXOrg <- crossPEDLinearReg(crossXOrg, c("xOrg"))
summary(lmXOrg)
## 
## Call:
## lm(formula = Q ~ pXOrg + pYOrg + pXNonOrg + pYNonOrg)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.65001 -0.23331 -0.06901  0.14514  1.17942 
## 
## Coefficients: (3 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -16.8432     0.5183  -32.49   <2e-16 ***
## pXOrg        -7.0398     0.1445  -48.73   <2e-16 ***
## pYOrg             NA         NA      NA       NA    
## pXNonOrg          NA         NA      NA       NA    
## pYNonOrg          NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3749 on 362 degrees of freedom
## Multiple R-squared:  0.8677, Adjusted R-squared:  0.8674 
## F-statistic:  2375 on 1 and 362 DF,  p-value: < 2.2e-16
lmXOrg <- crossPEDLinearReg(crossXOrg, c("yOrg"))
summary(lmXOrg)
## 
## Call:
## lm(formula = Q ~ pXOrg + pYOrg + pXNonOrg + pYNonOrg)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5019 -0.6540 -0.5033  0.5381  2.4173 
## 
## Coefficients: (3 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.5488     1.3545    7.05 9.14e-12 ***
## pXOrg             NA         NA      NA       NA    
## pYOrg         0.3226     0.3796    0.85    0.396    
## pXNonOrg          NA         NA      NA       NA    
## pYNonOrg          NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.03 on 362 degrees of freedom
## Multiple R-squared:  0.001991,   Adjusted R-squared:  -0.0007657 
## F-statistic: 0.7223 on 1 and 362 DF,  p-value: 0.396
lmXOrg <- crossPEDLinearReg(crossXOrg, c("xNonOrg"))
summary(lmXOrg)
## 
## Call:
## lm(formula = Q ~ pXOrg + pYOrg + pXNonOrg + pYNonOrg)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7368 -0.6850 -0.4753  0.6117  2.3207 
## 
## Coefficients: (3 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.9253     1.6707   3.547 0.000442 ***
## pXOrg             NA         NA      NA       NA    
## pYOrg             NA         NA      NA       NA    
## pXNonOrg     -0.6951     0.4693  -1.481 0.139438    
## pYNonOrg          NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.028 on 362 degrees of freedom
## Multiple R-squared:  0.006024,   Adjusted R-squared:  0.003278 
## F-statistic: 2.194 on 1 and 362 DF,  p-value: 0.1394
lmXOrg <- crossPEDLinearReg(crossXOrg, c("yNonOrg"))
summary(lmXOrg)
## 
## Call:
## lm(formula = Q ~ pXOrg + pYOrg + pXNonOrg + pYNonOrg)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5608 -0.6622 -0.4882  0.5950  2.4324 
## 
## Coefficients: (3 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.63069    1.35446   6.372 5.67e-10 ***
## pXOrg             NA         NA      NA       NA    
## pYOrg             NA         NA      NA       NA    
## pXNonOrg          NA         NA      NA       NA    
## pYNonOrg     0.06509    0.37955   0.171    0.864    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.031 on 362 degrees of freedom
## Multiple R-squared:  8.123e-05,  Adjusted R-squared:  -0.002681 
## F-statistic: 0.02941 on 1 and 362 DF,  p-value: 0.8639

Profit versus Price for the four products

XOrg.Qcoef = c(-16.8432,-7.0398)
YOrg.Qcoef = c(-15.7803,-6.7308)
XNonOrg.Qcoef = c(-15.8704,-6.7511)
YNonOrg.Qcoef = c(-15.4133,-6.657)

XOrg.Ccoef = c(-0.0006,0.8491)
YOrg.Ccoef = c(0.0003,0.8042)
XNonOrg.Ccoef = c(-0.0004,0.8434)
YNonOrg.Ccoef = c(0.0003,0.804)

profit <- function(price, Qcoef, Ccoef) {
  quantity <- exp(Qcoef[1]) * price^Qcoef[2]
  cost <- Ccoef[1] + Ccoef[2]*price
  return(quantity*(price-cost))
}
Xorg.Profit.func <- function(price) {
  return(profit(price,XOrg.Qcoef,XOrg.Ccoef))
}
Yorg.Profit.func <- function(price) {
  return(profit(price,YOrg.Qcoef,YOrg.Ccoef))
}
Xnonorg.Profit.func <- function(price) {
  return(profit(price,XNonOrg.Qcoef,XNonOrg.Ccoef))
}
Ynonorg.Profit.func <- function(price) {
  return(profit(price,YNonOrg.Qcoef,YNonOrg.Ccoef))
}
ggplot(data.frame(price = c(0.018, 0.035)), aes(price)) + 
  stat_function(fun = Xorg.Profit.func, aes(colour = "organic X")) + 
  stat_function(fun = Yorg.Profit.func, aes(colour = "organic Y")) + 
  stat_function(fun = Xnonorg.Profit.func, aes(colour = "inorganic X")) + 
  stat_function(fun = Ynonorg.Profit.func, aes(colour = "inorganic Y")) +
  scale_colour_manual("Lgend title", values = c("red", "blue", "green", "orange")) +
  labs(x = "price/oz", y = "profit/wk", title = "Profit per week versus price for each product")

Marginal Cost = Marginal Revenue

# Linear Model for C(q)
costQuantityLM <- function(data, CLASS, BRAND) {
  cost <- data %>% filter(class==CLASS, brand==BRAND) %>% select(c) %>% log() %>% as.matrix()
  quantity <- data %>% filter(class==CLASS, brand==BRAND) %>% select(oz) %>% log() %>% as.matrix()
  return(lm(cost~quantity))
}

# Linear Model for R(q)
revenueQuantityLM <- function(data, CLASS, BRAND) {
  revenue <- data %>% filter(class==CLASS, brand==BRAND) %>% select(p) %>% log() %>% as.matrix()
  quantity <- data %>% filter(class==CLASS, brand==BRAND) %>% select(oz) %>% log() %>% as.matrix()
  return(lm(revenue~quantity))
}

marginalValue <- function(point, coef) {
  return(exp(coef[1])*point^coef[2] + coef[2]*exp(coef[1])*point^(coef[2]-1)*point)
}

cQ_X_org <- costQuantityLM(data.ref, "organic", "X")
cQ_X_org$coefficients[1] <- cQ_X_org$coefficients[1]
cQ_X_org$coefficients[2] <- cQ_X_org$coefficients[2]/1.15
cQ_Y_org <- costQuantityLM(data.ref, "organic", "Y")
cQ_X_nonorg <- costQuantityLM(data.ref, "nonorganic", "X")
cQ_Y_nonorg <- costQuantityLM(data.ref, "nonorganic", "Y")

rQ_X_org <- revenueQuantityLM(data.ref, "organic", "X")
rQ_Y_org <- revenueQuantityLM(data.ref, "organic", "Y")
rQ_X_nonorg <- revenueQuantityLM(data.ref, "nonorganic", "X")
rQ_Y_nonorg <- revenueQuantityLM(data.ref, "nonorganic", "Y")

cQ_XO_Curve <- function(quantity) return(marginalValue(quantity,cQ_X_org$coefficients))
rQ_XO_Curve <- function(quantity) return(marginalValue(quantity,rQ_X_org$coefficients))

cQ_YO_Curve <- function(quantity) return(marginalValue(quantity,cQ_Y_org$coefficients))
rQ_YO_Curve <- function(quantity) return(marginalValue(quantity,rQ_Y_org$coefficients))

cQ_XI_Curve <- function(quantity) return(marginalValue(quantity,cQ_X_nonorg$coefficients))
rQ_XI_Curve <- function(quantity) return(marginalValue(quantity,rQ_X_nonorg$coefficients))

cQ_YI_Curve <- function(quantity) return(marginalValue(quantity,cQ_Y_nonorg$coefficients))
rQ_YI_Curve <- function(quantity) return(marginalValue(quantity,rQ_Y_nonorg$coefficients))

quantity.range = c(10, 100000)
ggplot(data.frame(quantity = quantity.range ), aes(quantity)) + 
  stat_function(fun = cQ_XO_Curve, aes(colour = "MC XO")) + 
  stat_function(fun = rQ_XO_Curve, aes(colour = "MR XO")) + 
  labs(y = "price/cost (per oz)", x="Quantity", title = "MC and MR for Organic X")

ggplot(data.frame(quantity = quantity.range ), aes(quantity)) + 
  stat_function(fun = cQ_YO_Curve, aes(colour = "MC YO")) + 
  stat_function(fun = rQ_YO_Curve, aes(colour = "MR YO")) + 
  labs(y = "price/cost (per oz)", x="Quantity", title = "MC and MR for Organic Y")

ggplot(data.frame(quantity = quantity.range ), aes(quantity)) + 
  stat_function(fun = cQ_XI_Curve, aes(colour = "MC XI")) + 
  stat_function(fun = rQ_XI_Curve, aes(colour = "MR XI")) + 
  labs(y = "price/cost (per oz)", x="Quantity", title = "MC and MR for NonOrganic X")

ggplot(data.frame(quantity = quantity.range ), aes(quantity)) + 
  stat_function(fun = cQ_YI_Curve, aes(colour = "MC YI")) + 
  stat_function(fun = rQ_YI_Curve, aes(colour = "MR YI")) + 
  labs(y = "price/cost (per oz)", x="Quantity", title = "MC and MR for NonOrganic Y")

MC = MR fake data

Reduced cost/oz response to an increase in quantity by 15% to reflect a greater fixed cost.

cQ_X_org <- costQuantityLM(data.ref, "organic", "X")
cQ_X_org$coefficients[1] <- cQ_X_org$coefficients[1]
cQ_X_org$coefficients[2] <- cQ_X_org$coefficients[2]/1.15
cQ_XO_Curve <- function(quantity) return(marginalValue(quantity,cQ_X_org$coefficients))
rQ_XO_Curve <- function(quantity) return(marginalValue(quantity,rQ_X_org$coefficients))
quantity.range = c(10, 100000)
ggplot(data.frame(quantity = quantity.range ), aes(quantity)) + 
  stat_function(fun = cQ_XO_Curve, aes(colour = "MC XO")) + 
  stat_function(fun = rQ_XO_Curve, aes(colour = "MR XO")) + 
  labs(y = "price/cost (per oz)", x="Quantity", title = "MC and MR for Organic X")

ggplot(data.frame(quantity = quantity.range ), aes(quantity)) + 
  stat_function(fun = cQ_XO_Curve, aes(colour = "MC XO")) + 
  stat_function(fun = rQ_XO_Curve, aes(colour = "MR XO")) + 
  labs(y = "price/cost (per oz)", x="Quantity", title = "MC and MR for Organic X") +
  scale_x_log10()

Multivariate cross-PED

crossPEDLM <- function(data.ref, BRAND="X", CLASS="organic") {
  quantity <- data.ref %>% filter(brand==BRAND, class==CLASS) %>% select(oz) %>% log() %>% as.matrix()
  hVal <- data.ref %>% filter(brand==BRAND, class==CLASS) %>% select(Hval) %>% log() %>% as.matrix()
  XO_org_price <- data.ref %>% filter(brand=="X", class=="organic") %>% select(p) %>% log() %>% as.matrix()
  YO_org_price <- data.ref %>% filter(brand=="Y", class=="organic") %>% select(p) %>% log() %>% as.matrix()
  XI_org_price <- data.ref %>% filter(brand=="X", class=="nonorganic") %>% select(p) %>% log() %>% as.matrix()
  YI_org_price <- data.ref %>% filter(brand=="Y", class=="nonorganic") %>% select(p) %>% log() %>% as.matrix()
  XO_org_deal <- data.ref %>% filter(brand=="X", class=="organic") %>% select(deal) %>% as.matrix()
  YO_org_deal <- data.ref %>% filter(brand=="Y", class=="organic") %>% select(deal) %>% as.matrix()
  XI_org_deal <- data.ref %>% filter(brand=="X", class=="nonorganic") %>% select(deal) %>% as.matrix()
  YI_org_deal <- data.ref %>% filter(brand=="Y", class=="nonorganic") %>% select(deal) %>% as.matrix()
  
  lm.crossPED <- lm(quantity~hVal + XO_org_price + YO_org_price + XI_org_price + YI_org_price + XO_org_deal + YO_org_deal + XI_org_deal + YI_org_deal)
}

Multivariate cross-PED for XO

crossPEDLMXO <- function(data.ref, BRAND="X", CLASS="organic") {
  quantity <- data.ref %>% filter(brand==BRAND, class==CLASS) %>% select(oz) %>% log() %>% as.matrix()
  hVal <- data.ref %>% filter(brand==BRAND, class==CLASS) %>% select(Hval) %>% log() %>% as.matrix()
  XO_org_price <- data.ref %>% filter(brand=="X", class=="organic") %>% select(p) %>% log() %>% as.matrix()
  YO_org_price <- data.ref %>% filter(brand=="Y", class=="organic") %>% select(p) %>% log() %>% as.matrix()
  XI_org_price <- data.ref %>% filter(brand=="X", class=="nonorganic") %>% select(p) %>% log() %>% as.matrix()
  YI_org_price <- data.ref %>% filter(brand=="Y", class=="nonorganic") %>% select(p) %>% log() %>% as.matrix()
  XO_org_deal <- data.ref %>% filter(brand=="X", class=="organic") %>% select(deal) %>% as.matrix()
  YO_org_deal <- data.ref %>% filter(brand=="Y", class=="organic") %>% select(deal) %>% as.matrix()
  XI_org_deal <- data.ref %>% filter(brand=="X", class=="nonorganic") %>% select(deal) %>% as.matrix()
  YI_org_deal <- data.ref %>% filter(brand=="Y", class=="nonorganic") %>% select(deal) %>% as.matrix()
  
  lm.crossPED <- lm(quantity~XO_org_price + XO_org_deal + YO_org_deal)
}
XO_cross_PED <- crossPEDLM(data.ref, "X", "organic")
summary(XO_cross_PED)
## 
## Call:
## lm(formula = quantity ~ hVal + XO_org_price + YO_org_price + 
##     XI_org_price + YI_org_price + XO_org_deal + YO_org_deal + 
##     XI_org_deal + YI_org_deal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35984 -0.11885 -0.03342  0.05692  1.16674 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.12343    1.31244  -0.094 0.925128    
## hVal          0.01133    0.02571   0.441 0.659719    
## XO_org_price -3.05819    0.19334 -15.818  < 2e-16 ***
## YO_org_price  0.66601    0.17098   3.895 0.000117 ***
## XI_org_price -0.12149    0.18959  -0.641 0.522066    
## YI_org_price  0.23843    0.17355   1.374 0.170350    
## XO_org_deal   1.43078    0.06176  23.166  < 2e-16 ***
## YO_org_deal   0.02100    0.06157   0.341 0.733207    
## XI_org_deal  -0.04871    0.06125  -0.795 0.426987    
## YI_org_deal   0.08843    0.06172   1.433 0.152804    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2322 on 354 degrees of freedom
## Multiple R-squared:  0.9504, Adjusted R-squared:  0.9491 
## F-statistic: 753.3 on 9 and 354 DF,  p-value: < 2.2e-16
XO_cross_PED <- crossPEDLMXO(data.ref, "X", "organic")
summary(XO_cross_PED)
## 
## Call:
## lm(formula = quantity ~ XO_org_price + XO_org_deal + YO_org_deal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34544 -0.13032 -0.03795  0.05013  1.18062 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3.16105    0.68321  -4.627 5.19e-06 ***
## XO_org_price -3.13756    0.19424 -16.153  < 2e-16 ***
## XO_org_deal   1.39063    0.06125  22.704  < 2e-16 ***
## YO_org_deal  -0.18046    0.03159  -5.712 2.34e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2362 on 360 degrees of freedom
## Multiple R-squared:  0.9478, Adjusted R-squared:  0.9473 
## F-statistic:  2178 on 3 and 360 DF,  p-value: < 2.2e-16

Multivariate cross-PED for YO

crossPEDLMYO <- function(data.ref, BRAND="Y", CLASS="organic") {
  quantity <- data.ref %>% filter(brand==BRAND, class==CLASS) %>% select(oz) %>% log() %>% as.matrix()
  hVal <- data.ref %>% filter(brand==BRAND, class==CLASS) %>% select(Hval) %>% log() %>% as.matrix()
  XO_org_price <- data.ref %>% filter(brand=="X", class=="organic") %>% select(p) %>% log() %>% as.matrix()
  YO_org_price <- data.ref %>% filter(brand=="Y", class=="organic") %>% select(p) %>% log() %>% as.matrix()
  XI_org_price <- data.ref %>% filter(brand=="X", class=="nonorganic") %>% select(p) %>% log() %>% as.matrix()
  YI_org_price <- data.ref %>% filter(brand=="Y", class=="nonorganic") %>% select(p) %>% log() %>% as.matrix()
  XO_org_deal <- data.ref %>% filter(brand=="X", class=="organic") %>% select(deal) %>% as.matrix()
  YO_org_deal <- data.ref %>% filter(brand=="Y", class=="organic") %>% select(deal) %>% as.matrix()
  XI_org_deal <- data.ref %>% filter(brand=="X", class=="nonorganic") %>% select(deal) %>% as.matrix()
  YI_org_deal <- data.ref %>% filter(brand=="Y", class=="nonorganic") %>% select(deal) %>% as.matrix()
  
  lm.crossPED <- lm(quantity~XO_org_price + YO_org_price + YO_org_deal)
}
YO_cross_PED <- crossPEDLM(data.ref, "Y", "organic")
summary(YO_cross_PED)
## 
## Call:
## lm(formula = quantity ~ hVal + XO_org_price + YO_org_price + 
##     XI_org_price + YI_org_price + XO_org_deal + YO_org_deal + 
##     XI_org_deal + YI_org_deal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.94464 -0.64035  0.00247  0.59479  2.31866 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.366189   5.332936   0.256   0.7980    
## hVal         -0.079762   0.104471  -0.763   0.4457    
## XO_org_price  0.557819   0.785592   0.710   0.4781    
## YO_org_price -3.954382   0.694750  -5.692 2.64e-08 ***
## XI_org_price  1.372694   0.770372   1.782   0.0756 .  
## YI_org_price  0.215070   0.705186   0.305   0.7606    
## XO_org_deal  -0.092792   0.250961  -0.370   0.7118    
## YO_org_deal   1.176298   0.250179   4.702 3.70e-06 ***
## XI_org_deal   0.455553   0.248870   1.830   0.0680 .  
## YI_org_deal   0.005135   0.250776   0.020   0.9837    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9436 on 354 degrees of freedom
## Multiple R-squared:  0.535,  Adjusted R-squared:  0.5232 
## F-statistic: 45.26 on 9 and 354 DF,  p-value: < 2.2e-16
YO_cross_PED <- crossPEDLMYO(data.ref, "Y", "organic")
summary(YO_cross_PED)
## 
## Call:
## lm(formula = quantity ~ XO_org_price + YO_org_price + YO_org_deal)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8533 -0.6560  0.0318  0.5870  2.7387 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -3.4693     2.6189  -1.325   0.1861    
## XO_org_price   0.7786     0.3634   2.143   0.0328 *  
## YO_org_price  -3.9990     0.6721  -5.950 6.36e-09 ***
## YO_org_deal    1.1654     0.2438   4.780 2.56e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9412 on 360 degrees of freedom
## Multiple R-squared:  0.5296, Adjusted R-squared:  0.5256 
## F-statistic: 135.1 on 3 and 360 DF,  p-value: < 2.2e-16

Multivariate cross-PED for XI

crossPEDLMXI <- function(data.ref, BRAND="X", CLASS="nonorganic") {
  quantity <- data.ref %>% filter(brand==BRAND, class==CLASS) %>% select(oz) %>% log() %>% as.matrix()
  hVal <- data.ref %>% filter(brand==BRAND, class==CLASS) %>% select(Hval) %>% log() %>% as.matrix()
  XO_org_price <- data.ref %>% filter(brand=="X", class=="organic") %>% select(p) %>% log() %>% as.matrix()
  YO_org_price <- data.ref %>% filter(brand=="Y", class=="organic") %>% select(p) %>% log() %>% as.matrix()
  XI_org_price <- data.ref %>% filter(brand=="X", class=="nonorganic") %>% select(p) %>% log() %>% as.matrix()
  YI_org_price <- data.ref %>% filter(brand=="Y", class=="nonorganic") %>% select(p) %>% log() %>% as.matrix()
  XO_org_deal <- data.ref %>% filter(brand=="X", class=="organic") %>% select(deal) %>% as.matrix()
  YO_org_deal <- data.ref %>% filter(brand=="Y", class=="organic") %>% select(deal) %>% as.matrix()
  XI_org_deal <- data.ref %>% filter(brand=="X", class=="nonorganic") %>% select(deal) %>% as.matrix()
  YI_org_deal <- data.ref %>% filter(brand=="Y", class=="nonorganic") %>% select(deal) %>% as.matrix()
  
  lm.crossPED <- lm(quantity~XI_org_price + YI_org_price + XI_org_deal)
}
XI_cross_PED <- crossPEDLM(data.ref, "X", "nonorganic")
summary(XI_cross_PED)
## 
## Call:
## lm(formula = quantity ~ hVal + XO_org_price + YO_org_price + 
##     XI_org_price + YI_org_price + XO_org_deal + YO_org_deal + 
##     XI_org_deal + YI_org_deal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.40513 -0.09481 -0.01748  0.05342  1.27001 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.80247    1.14292  -0.702  0.48307    
## hVal          0.01148    0.02239   0.513  0.60838    
## XO_org_price -0.05371    0.16836  -0.319  0.74991    
## YO_org_price  0.15148    0.14889   1.017  0.30968    
## XI_org_price -3.01591    0.16510 -18.267  < 2e-16 ***
## YI_org_price  0.45680    0.15113   3.023  0.00269 ** 
## XO_org_deal  -0.03101    0.05378  -0.577  0.56463    
## YO_org_deal   0.05907    0.05362   1.102  0.27136    
## XI_org_deal   1.45074    0.05334  27.200  < 2e-16 ***
## YI_org_deal  -0.02728    0.05374  -0.508  0.61209    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2022 on 354 degrees of freedom
## Multiple R-squared:  0.9453, Adjusted R-squared:  0.9439 
## F-statistic: 679.8 on 9 and 354 DF,  p-value: < 2.2e-16
XI_cross_PED <- crossPEDLMXI(data.ref, "X", "nonorganic")
summary(XI_cross_PED)
## 
## Call:
## lm(formula = quantity ~ XI_org_price + YI_org_price + XI_org_deal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38903 -0.09652 -0.02595  0.05932  1.27352 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.94038    0.65699  -1.431    0.153    
## XI_org_price -3.02961    0.16307 -18.579  < 2e-16 ***
## YI_org_price  0.53499    0.07467   7.165 4.44e-12 ***
## XI_org_deal   1.44573    0.05236  27.612  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2013 on 360 degrees of freedom
## Multiple R-squared:  0.9449, Adjusted R-squared:  0.9444 
## F-statistic:  2057 on 3 and 360 DF,  p-value: < 2.2e-16

Multivariate cross-PED for YI

crossPEDLMYI <- function(data.ref, BRAND="Y", CLASS="nonorganic") {
  quantity <- data.ref %>% filter(brand==BRAND, class==CLASS) %>% select(oz) %>% log() %>% as.matrix()
  hVal <- data.ref %>% filter(brand==BRAND, class==CLASS) %>% select(Hval) %>% log() %>% as.matrix()
  XO_org_price <- data.ref %>% filter(brand=="X", class=="organic") %>% select(p) %>% log() %>% as.matrix()
  YO_org_price <- data.ref %>% filter(brand=="Y", class=="organic") %>% select(p) %>% log() %>% as.matrix()
  XI_org_price <- data.ref %>% filter(brand=="X", class=="nonorganic") %>% select(p) %>% log() %>% as.matrix()
  YI_org_price <- data.ref %>% filter(brand=="Y", class=="nonorganic") %>% select(p) %>% log() %>% as.matrix()
  XO_org_deal <- data.ref %>% filter(brand=="X", class=="organic") %>% select(deal) %>% as.matrix()
  YO_org_deal <- data.ref %>% filter(brand=="Y", class=="organic") %>% select(deal) %>% as.matrix()
  XI_org_deal <- data.ref %>% filter(brand=="X", class=="nonorganic") %>% select(deal) %>% as.matrix()
  YI_org_deal <- data.ref %>% filter(brand=="Y", class=="nonorganic") %>% select(deal) %>% as.matrix()
  
  lm.crossPED <- lm(quantity~XI_org_price + YI_org_price + YI_org_deal)
}
YI_cross_PED <- crossPEDLM(data.ref, "Y", "nonorganic")
summary(YI_cross_PED)
## 
## Call:
## lm(formula = quantity ~ hVal + XO_org_price + YO_org_price + 
##     XI_org_price + YI_org_price + XO_org_deal + YO_org_deal + 
##     XI_org_deal + YI_org_deal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.57701 -0.56374  0.00957  0.54617  2.74225 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.3969129  5.1438843  -0.272    0.786    
## hVal         -0.1313873  0.1007674  -1.304    0.193    
## XO_org_price -0.5463521  0.7577431  -0.721    0.471    
## YO_org_price -0.1401418  0.6701215  -0.209    0.834    
## XI_org_price  1.0388633  0.7430626   1.398    0.163    
## YI_org_price -2.9248065  0.6801875  -4.300 2.21e-05 ***
## XO_org_deal  -0.0610505  0.2420646  -0.252    0.801    
## YO_org_deal   0.0006412  0.2413103   0.003    0.998    
## XI_org_deal   0.1519458  0.2400471   0.633    0.527    
## YI_org_deal   1.5008699  0.2418856   6.205 1.53e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9101 on 354 degrees of freedom
## Multiple R-squared:  0.5554, Adjusted R-squared:  0.5441 
## F-statistic: 49.13 on 9 and 354 DF,  p-value: < 2.2e-16
YI_cross_PED <- crossPEDLMYI(data.ref, "Y", "nonorganic")
summary(YI_cross_PED)
## 
## Call:
## lm(formula = quantity ~ XI_org_price + YI_org_price + YI_org_deal)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4585 -0.5614  0.0186  0.5795  2.8311 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.2907     2.7120  -0.107    0.915    
## XI_org_price   0.6184     0.4151   1.490    0.137    
## YI_org_price  -2.9475     0.6710  -4.393 1.47e-05 ***
## YI_org_deal    1.5197     0.2387   6.366 5.92e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9072 on 360 degrees of freedom
## Multiple R-squared:  0.5507, Adjusted R-squared:  0.547 
## F-statistic: 147.1 on 3 and 360 DF,  p-value: < 2.2e-16
summary_XO <- data.ref %>% filter(brand=="X", class=="organic") %>% select(oz, p, c, deal)
summary_XO <- mutate(summary_XO, profit=oz*(p-c))
summary_XO$deal[summary_XO$deal == 0] <- "Without"
summary_XO$deal[summary_XO$deal == 1] <- "With"
ggplot(summary_XO, aes(x=deal, y=profit)) + geom_boxplot() + labs(title="XO profit with and without a deal", x="Deal", y = "Profit") 

ggplot(summary_XO, aes(x=deal, y=oz)) + geom_boxplot() + labs(title="XO profit with and without a deal", x="Deal", y = "Oz") 

summary_YO <- data.ref %>% filter(brand=="Y", class=="organic") %>% select(oz, p, c, deal)
summary_YO <- mutate(summary_YO, profit=oz*(p-c))
summary_YO$deal[summary_YO$deal == 0] <- "Without"
summary_YO$deal[summary_YO$deal == 1] <- "With"
ggplot(summary_YO, aes(x=deal, y=profit)) + geom_boxplot() + labs(title="XO profit with and without a deal", x="Deal", y = "Profit") 

ggplot(summary_YO, aes(x=deal, y=oz)) + geom_boxplot() + labs(title="XO profit with and without a deal", x="Deal", y = "Oz") 

summary_XI <- data.ref %>% filter(brand=="X", class=="nonorganic") %>% select(oz, p, c, deal)
summary_XI <- mutate(summary_XI, profit=oz*(p-c))
summary_XI$deal[summary_XI$deal == 0] <- "Without"
summary_XI$deal[summary_XI$deal == 1] <- "With"
ggplot(summary_XI, aes(x=deal, y=profit)) + geom_boxplot() + labs(title="XO profit with and without a deal", x="Deal", y = "Profit") 

ggplot(summary_XI, aes(x=deal, y=oz)) + geom_boxplot() + labs(title="XO profit with and without a deal", x="Deal", y = "Oz") 

summary_YI <- data.ref %>% filter(brand=="Y", class=="nonorganic") %>% select(oz, p, c, deal)
summary_YI <- mutate(summary_YI, profit=oz*(p-c))
summary_YI$deal[summary_YI$deal == 0] <- "Without"
summary_YI$deal[summary_YI$deal == 1] <- "With"
ggplot(summary_YI, aes(x=deal, y=profit)) + geom_boxplot() + labs(title="XO profit with and without a deal", x="Deal", y = "Profit") 

ggplot(summary_YI, aes(x=deal, y=oz)) + geom_boxplot() + labs(title="XO profit with and without a deal", x="Deal", y = "Oz")